Instrucciones generales

Genere un R Markdown en donde se muestre el código utilizado y la salida de este. Utilice los paquetes asociados al tidyverse para manipular los datos (tidyr, dplyr, …) y generar las visualizaciones (ggplot2), así como el paquete lubridate para manipular fechas.

Leer los datos

  1. Cargue los paquetes que vaya a utilizar en el análisis.

    library(tidyverse)
    ## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
    ## v ggplot2 3.3.3     v purrr   0.3.4
## v tibble  3.1.1     v dplyr   1.0.5
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.1
    ## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
    library(lubridate)
    ## 
## Attaching package: 'lubridate'
    ## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
    library(plotly)
    ## 
## Attaching package: 'plotly'
    ## The following object is masked from 'package:ggplot2':
## 
##     last_plot
    ## The following object is masked from 'package:stats':
## 
##     filter
    ## The following object is masked from 'package:graphics':
## 
##     layout
    library(DT)

  2. Lea los datos del archivo qd.csv, asígnelos a un objeto y escoja dos de las series de caudales promedio diarios que se presentan en las columnas de la tabla de datos (denotadas eh1, eh2, …). Sugerencia: utilice el paquete readr.

    # inserte acá su código
#setwd("C:/Users/Jesus/Dropbox/UCR/2021/Computacional/Lab 2")
setwd("C:/Users/rzamoram/OneDrive - Intel Corporation/Documents/Machine Learning/Estadistica Computacional")

qd <- readr::read_csv("qd.csv",col_types = cols(fecha = col_date(format = "%Y-%m-%d")))

Serie de caudales diarios

  1. Calcule el número de datos faltantes de cada serie.

    # inserte acá su código

  apply(qd, 2, function(x) sum(is.na(x)))
    ## fecha   eh1   eh2   eh3   eh4   eh5   eh6   eh7   eh8   eh9  eh10  eh11  eh12 
##     0  1887  1947  2037  1509  2544  1689  2402  2084  2206  1654  1314  2687 
##  anno 
##     0

  2. Visualice las series de caudales diarios mediante gráficos de línea (i.e., dos gráficos independientes). Recuerde nombrar ejes, título del gráfico, etc.

    Para las siguientes preguntas utilice únicamente 3 años de datos consecutivos de cada estación hidrológica

  3. Visualice las series de caudales diarios mediante un único gráfico pero diferenciando cada una por color distinto.

    # inserte acá su código
  qd1 <- qd %>% 
  filter(year(fecha) >2016) %>%
  pivot_longer(cols = starts_with("eh"),names_to = "caudal",names_prefix = "eh",values_drop_na = TRUE) %>%
  mutate(caudal = as.factor(as.numeric(caudal))) %>% 
  ggplot(aes( x = fecha, y = value,color = caudal)) +
  geom_line()  + 
  labs(
    x = "Fecha",
    y = "Volumen diario del caudal",
    colour = "Caudal",
    title = "Series de volumen del caudal diario de los últimos 3 años") +
  scale_colour_brewer(type = "seq", palette = "Paired")+
  theme_bw() + 
  theme(
plot.title = element_text(face = "bold", size = 12),
legend.background = element_rect(fill = "grey70", size = 4, colour = "white"),
legend.justification = c(0, -1),
legend.position = c(2, -1),
axis.ticks = element_line(colour = "black", size = 0.2),
panel.grid.major = element_line(colour = "black", size = 0.2),
panel.grid.minor = element_blank()) 
  ggplotly(qd1)

  4. Visualice las series de caudales diarios mediante páneles independientes y alineados (un único gráfico). Sugerencia 1: use alguna función de tipo facet_*(). Sugerencia 2: de ser necesario, utilice el argumento scales de las funciones facet_*().

    # inserte acá su código
  qd1 <- qd %>% 
  filter(year(fecha) >2016) %>%
  pivot_longer(cols = starts_with("eh"),names_to = "caudal",names_prefix = "eh",values_drop_na = TRUE) %>%  #NO VER NAS
  mutate(caudal = as.numeric(caudal)) %>% 
  ggplot(aes( x = fecha, y = value)) +
  geom_line(color = "orange")  + 
  labs(
    x = "Fecha",
    y = "Volumen diario del caudal",
    colour = "Caudal",
    title = "Series de volumen del caudal diario de los últimos 3 años") +
  scale_colour_brewer(type = "seq", palette = "Paired")+
  theme_bw() + 
  theme(
plot.title = element_text(face = "bold", size = 12),
legend.background = element_rect(fill = "grey70", size = 4, colour = "white"),
legend.justification = c(0, -1),
legend.position = c(2, -1),
axis.ticks = element_line(colour = "black", size = 0.2),
panel.grid.major = element_line(colour = "black", size = 0.2),
panel.grid.minor = element_blank()) +
facet_wrap(caudal~.,ncol=3,scales = "free") 
  ggplotly(qd1)

  5. Visualice los periodos de datos faltantes de la serie. Puede ser ya sea generando bandas verticales traslucidas en los sitios donde se tienen faltantes de datos, graficando líneas verticales en dichas puntos o utilizando algún código de colores para la líinea en dichos tramos. Construya un arreglo en el cual se tengan las estaciones en las columnas y los años seleccionados en las filas.

  • Construya un arreglo en el cual se tengan las estaciones en las columnas y los años seleccionados en las filas.

    # inserte acá su código
  qd1 <- qd %>% 
  filter(year(fecha) >2016) %>%
  pivot_longer(cols = starts_with("eh"),names_to = "caudal",names_prefix = "eh",values_drop_na = FALSE) %>% #presentar NAs
  mutate(caudal = as.numeric(caudal),caudal=as.factor(caudal)) 
  nans <- qd1 %>% filter(is.na(value)==TRUE) %>% mutate(caudal=as.factor(caudal))

  qd1 %>% filter(is.na(value)==TRUE) %>%
  mutate(caudal=as.factor(caudal)) %>% 
  group_by(anno, caudal) %>% 
  summarise(Nans =n()) %>% 
  pivot_wider(names_from = caudal, values_from = Nans, values_fill =0 ) %>% 
  datatable()
    ## `summarise()` has grouped output by 'anno'. You can override using the `.groups` argument.

  • Visualización

    # inserte acá su código
  qd1 <- qd %>% 
  filter(year(fecha) >2016) %>%
  pivot_longer(cols = starts_with("eh"),names_to = "caudal",names_prefix = "eh",values_drop_na = FALSE) %>%
  mutate(caudal = as.numeric(caudal),caudal=as.factor(caudal)) 
  nans <- qd1 %>% filter(is.na(value)==TRUE) %>% mutate(caudal=as.factor(caudal))

  qd1 %>% 
  ggplot(aes( x = fecha, y = value, group = caudal)) +
  geom_line(color = "red")  + 
  labs(
    x = "Fecha",
    y = "Volumen diario del caudal",
    colour = "Caudal",
    title = "Visualización de los períodos de datos faltantes") +
  geom_vline(aes(xintercept=fecha), data=nans,color = 'black',alpha=0.05) +
  theme_bw() + 
  theme(
plot.title = element_text(face = "bold", size = 12),
legend.background = element_rect(fill = "grey70", size = 4, colour = "white"),
legend.justification = c(0, -1),
legend.position = c(2, -1),
axis.ticks = element_line(colour = "grey70", size = 0.2),
panel.grid.major = element_line(colour = "grey70", size = 0.2),
panel.grid.minor = element_blank()) +
facet_wrap(caudal~.,ncol=3,scales = "free")
    ## Warning: Removed 183 row(s) containing missing values (geom_path).

    Para las siguientes preguntas utilice todo el registro de las estaciones hidrológicas

  1. Visualice la función de distribución empírica de cada estación. Sugerencia: utilice una función del paquete ggplot2 que termina con _ecdf().

    # inserte acá su código
  qd1 <- qd %>% 
  filter(year(fecha) >2016) %>%
  pivot_longer(cols = starts_with("eh"),names_to = "caudal",names_prefix = "eh",values_drop_na = TRUE) %>%
  mutate(caudal = as.numeric(caudal), caudal = as.factor(caudal)) %>% 
  ggplot(aes( x = fecha, y = value, color = caudal)) +
  stat_ecdf() +
  labs(
    x = "Fecha",
    y = "Función distribución empírica de cada Estación", 
    colour = "Caudal",
    title = "Función distribución empírica de cada Estación desde 2017 a 2019") +
  theme_bw() + 
  theme(
plot.title = element_text(face = "bold", size = 12),
legend.background = element_rect(fill = "grey70", size = 4, colour = "white"),
legend.justification = c(0, -1),
legend.position = c(2, -1),
axis.ticks = element_line(colour = "grey70", size = 0.2),
panel.grid.major = element_line(colour = "grey70", size = 0.2),
panel.grid.minor = element_blank()) +
facet_wrap(caudal~.,ncol=3,scales = "free") 
ggplotly(qd1)

  2. Para los registros comunes, realice un gráfico de puntos (scatterplot) en donde los ejes representen los caudales promedio diarios de cada una de las estaciones.

    # inserte acá su código
  qd %>% 
  filter(year(fecha) >2016) %>%
  pivot_longer(cols = starts_with("eh"),names_to = "caudal",names_prefix = "eh",values_drop_na = TRUE) %>%
  mutate(caudal = as.numeric(caudal), caudal = as.factor(caudal)) %>% 
  group_by(caudal) %>% 
  summarise(value = mean(value)) %>% 
  ggplot(aes( x = caudal, y = value)) +
  geom_point(color = "red", size = 3) +
  theme_bw() + 
  theme(
plot.title = element_text(face = "bold", size = 12),
legend.background = element_rect(fill = "grey70", size = 4, colour = "white"),
legend.justification = c(0, -1),
legend.position = c(2, -1),
axis.ticks = element_line(colour = "black", size = 0.2),
panel.grid.major = element_line(colour = "black", size = 0.2),
panel.grid.minor = element_blank())

Serie de caudales a nivel mensual

  1. Calcule el número de datos faltantes de cada serie por cada mes y año. Sugerencia: cree una nueva columna con una indicadora que identifique si el dato está ausente o no utilizando is.na().

    # Generar la columna de mes para cada observación
base100 <- qd   
qd$mes <-(month(qd$fecha))

base2 <- qd %>% pivot_longer(c(eh4,eh11), names_to = "series", values_to = "valor")

#Contar los NAs
base2$na <- rep(NA, nrow(base2))

  for (i in 1:nrow(base2)) {
if (is.na(base2$valor[i]) == TRUE) {
base2$na[i] <- 1
 } else
base2$na[i] <- 0
  } 

#Calcular los NAs por mes
tabla1 <- base2 %>%  group_by(series, mes) %>% summarise(Faltantes = sum(na)) %>% ungroup()
    ## `summarise()` has grouped output by 'series'. You can override using the `.groups` argument.
    datatable(tabla1)
    #Ahora por año

tabla2 <- base2 %>%  group_by(series, anno) %>% summarise(Faltantes = sum(na)) %>% ungroup()
    ## `summarise()` has grouped output by 'series'. You can override using the `.groups` argument.
    datatable(tabla2)

  2. Visualice el número de datos faltantes como una serie de datos. ¿Existe algún patrón en el cual se concentren la mayor cantidad de datos faltantes?

    gg1<-tabla1 %>% ggplot() + geom_line(aes(x=factor(mes), 
                       y=Faltantes,
                       col=series, group = 1))+
labs(title="Datos faltantes por mes de las Series de caudales eh4 y eh11", 
     x="mes",
     y="Datos faltantes",
     col= "Serie")+
  theme_bw()+
  facet_grid(~series)+
  theme(
plot.title = element_text(face = "bold", size = 12),
legend.background = element_rect(fill = "grey70", size = 4, colour = "white"),
legend.justification = c(0, -1),
legend.position = c(2, -1),
axis.ticks = element_line(colour = "grey70", size = 0.2),
panel.grid.major = element_line(colour = "grey70", size = 0.2),
panel.grid.minor = element_blank()) 

ggplotly(gg1)

  3. Calcule la serie de caudales promedio mensuales e incluya el número de datos con el cual se está calculando cada promedio.

    base3 <- drop_na(base2)  

base4 <-  base3 %>% group_by(mes, series) %>% summarise(promedio= mean(valor, na.rm=TRUE), n =n())  
    ## `summarise()` has grouped output by 'mes'. You can override using the `.groups` argument.
    base4$mes <-month(base4$mes)
base4 <- arrange(base4, series)
datatable(base4)

  4. Visualice las series de caudales mensuales mediante páneles independientes (una única serie por pánel).

    gg2<-base4 %>% ggplot() + geom_line(aes(x=mes, 
                       y=promedio,
                       col=series, group=1))+
labs(title="Caudales mensuales", 
     x="Mes",
     y="Caudal",
     col= "Serie")+ facet_grid(~series)+
  theme_bw()+
  theme(
plot.title = element_text(face = "bold", size = 12),
legend.background = element_rect(fill = "grey70", size = 4, colour = "white"),
legend.justification = c(0, -1),
legend.position = c(2, -1),
axis.ticks = element_line(colour = "grey70", size = 0.2),
panel.grid.major = element_line(colour = "grey70", size = 0.2),
panel.grid.minor = element_blank()) 

ggplotly(gg2)

  5. Cree una tabla por estación en la cual las columnas representen los meses y las filas los años. El valor a desplegar en cada celda sería el caudal promedio mensual. Sugerencia: utilice pivot_wider() o pivot_longer().

    qd$mes <-(month(ymd(qd$fecha), label = TRUE, abbr = FALSE))

tabla3 <- qd %>% pivot_wider(id_cols=anno, names_from = mes, values_from = eh4, values_fn =mean)
datatable(tabla3)
    tabla4 <- qd %>% pivot_wider(id_cols=anno, names_from = mes, values_from = eh11, values_fn =mean)
datatable(tabla4)

  6. Visualice las series de caudales mensuales.

    qd$mes <- as.character(qd$mes) 

# Para EH4    
mes <- qd %>% group_by(mes, anno) %>% summarise(v1=mean(eh4)) %>% ungroup()
    ## `summarise()` has grouped output by 'mes'. You can override using the `.groups` argument.
      gg3<-mes %>% ggplot() + geom_line(aes(x=anno, 
                       y=v1, col=mes))+
labs(title="Promedio mensual de EH4", 
     x="Fecha",
     y="Caudal diario",
     col= "Serie")+ theme_bw()+
theme(
plot.title = element_text(face = "bold", size = 12),
legend.background = element_rect(fill = "grey70", size = 4, colour = "white"),
legend.justification = c(0, -1),
legend.position = c(2, -1),
axis.ticks = element_line(colour = "grey70", size = 0.2),
panel.grid.major = element_line(colour = "grey70", size = 0.2),
panel.grid.minor = element_blank()) 

  ggplotly(gg3) 
      mes2 <- qd %>% group_by(mes, anno) %>% summarise(v1=mean(eh11)) %>% ungroup()
    ## `summarise()` has grouped output by 'mes'. You can override using the `.groups` argument.
      gg4<-mes2 %>% ggplot() + geom_line(aes(x=anno, 
                       y=v1, col=mes))+
labs(title="Serie promedio mensual eh11", 
     x="Fecha",
     y="Caudal diario",
     col= "Serie")+theme_bw()+
theme(
plot.title = element_text(face = "bold", size = 12),
legend.background = element_rect(fill = "grey70", size = 4, colour = "white"),
legend.justification = c(0, -1),
legend.position = c(2, -1),
axis.ticks = element_line(colour = "grey70", size = 0.2),
panel.grid.major = element_line(colour = "grey70", size = 0.2),
panel.grid.minor = element_blank()) 

  ggplotly(gg4)

  7. Visualice el régimen estacional, y su variabilidad, mediante un boxplot o un violin plot.

    gg5<-base2 %>% ggplot() + geom_boxplot(aes(x=series, y=valor)) +
  labs(title="Regimen Estacional", x="Estación", y="Caudal")+theme_bw()+theme(
plot.title = element_text(face = "bold", size = 12),
legend.background = element_rect(fill = "grey70", size = 4, colour = "white"),
legend.justification = c(0, -1),
legend.position = c(2, -1),
axis.ticks = element_line(colour = "grey70", size = 0.2),
panel.grid.major = element_line(colour = "grey70", size = 0.2),
panel.grid.minor = element_blank()) 

ggplotly(gg5)      
    ## Warning: Removed 2823 rows containing non-finite values (stat_boxplot).

Análisis varios

  1. Calcule la serie de caudales promedio anuales.
base13 <- drop_na(base2)  
    
base14 <-  base13 %>% group_by(anno, series) %>% summarise(promedio2= mean(valor, na.rm=TRUE), n =n())  
## `summarise()` has grouped output by 'anno'. You can override using the `.groups` argument.
#base14$anno <-year(base14$anno)
base14 <- arrange(base14, series)
datatable(base14)

  1. Calcule la serie de anual de caudales máximos diarios

    base23 <- drop_na(base2)  

base24 <-  base23 %>% group_by(anno, series) %>% summarise(maximo= max(valor, na.rm=TRUE), n =n())  
    ## `summarise()` has grouped output by 'anno'. You can override using the `.groups` argument.
    #base14$anno <-year(base14$anno)
base24 <- arrange(base24, series)
datatable(base24)

  2. Visualice la serie de caudales promedio diarios junto con los máximos anuales.

promax <- ggplot() + 
  geom_line(data = base14, aes(x = anno, y = promedio2, col=series, group=1)) +
  geom_line(data = base24, aes(x = anno, y = maximo, col=series, group=1)) +
    labs(title="Promedio de caudales diarios anuales junto con maximos anuales", 
         x="Fecha",
         y="Caudal diario",
         col= "Serie")+ theme_bw()+
    theme(
    plot.title = element_text(face = "bold", size = 12),
    legend.background = element_rect(fill = "grey70", size = 4, colour = "white"),
    legend.justification = c(0, -1),
    legend.position = c(2, -1),
    axis.ticks = element_line(colour = "grey70", size = 0.2),
    panel.grid.major = element_line(colour = "grey70", size = 0.2),
    panel.grid.minor = element_blank()) 

    

  ggplotly(promax)

  1. Con base en la función de distribución empírica, estime cuál es el caudal que se excede únicamente 5% del tiempo y extraiga los caudales que se encuentran por encima de este umbral.

    #qd5<-qd[,2:13]
qd7<-base100[base100$eh4 > quantile(base14$promedio2,prob=0.95, na.rm = TRUE),]
head(qd7[c("fecha", "eh4")])
    ## # A tibble: 6 x 2
##   fecha        eh4
##   <date>     <dbl>
## 1 1965-01-06  1.11
## 2 1965-01-07  2.27
## 3 1965-01-08  1.46
## 4 1965-01-18  4.50
## 5 1965-01-19  4.18
## 6 1965-01-20  2.39
    qd8<-base100[base100$eh11 > quantile(base14$promedio2,prob=0.95, na.rm = TRUE),]
head(qd8[c("fecha", "eh11")])
    ## # A tibble: 6 x 2
##   fecha       eh11
##   <date>     <dbl>
## 1 1965-06-24 0.952
## 2 1965-06-25 1.13 
## 3 1965-09-18 0.935
## 4 1965-09-20 1.31 
## 5 1965-09-21 1.93 
## 6 1965-09-22 1.33

  2. Visualice la serie de caudales promedio diarios junto con los valores que exceden el umbral estimado en el punto anterior.

promax <- ggplot() + 
  geom_line(data = base14, aes(x = anno, y = promedio2, col=series, group=1)) +
  geom_point(data = qd7, aes(x = anno, y = eh4, group=1, na.rm=T)) +
  geom_point(data = qd8, aes(x = anno, y = eh11, group=1, na.rm=T)) +
    labs(title="caudales promedio diarios junto con los valores que exceden el umbral", 
         x="Fecha",
         y="Caudal diario",
         col= "Serie")+ theme_bw()+
    theme(
    plot.title = element_text(face = "bold", size = 12),
    legend.background = element_rect(fill = "grey70", size = 4, colour = "white"),
    legend.justification = c(0, -1),
    legend.position = c(2, -1),
    axis.ticks = element_line(colour = "grey70", size = 0.2),
    panel.grid.major = element_line(colour = "grey70", size = 0.2),
    panel.grid.minor = element_blank()) 
## Warning: Ignoring unknown aesthetics: na.rm

## Warning: Ignoring unknown aesthetics: na.rm
  ggplotly(promax)